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1. Introduction 


The purpose of this paper is a calculation of the stationary space energy 
distribution of slow neutrons in a moderator at all distances from the source. 

Exact methods of calculating the density of slowed-down neutrons have been 
made available by the derivation of the general equation governing the diffu- 
sion and slowing-down of neutrons?, i.e. the transport equation, and its solu- 
tion for the case of an isotropic point source in an infinite medium by means 
of a Fourier-Laplace transformation?. The elementary age-theory? gives a general 
formula for the slow neutron density, which is assumed to be valid at small dis- 
tances r from the source, for a slowly varying mean free path A(F) of the 
neutrons (H = neutron energy), and best for large values of M (= mass of the 
scattering atomic nuclei measured in units of the neutron mass). In other 
cases, however, and always for large r, this formula may give considerable 
errors. Therefore it is important to make a calculation of the neutron density 
which is based on the general theory, and for this purpose the Fourier-Laplace 
transform of the transport equation seems to be the most convenient start- 
ing-point. 

So far, most calculations of the distribution of slowed-down neutrons out- 
side the range of validity of age-theory have been made assuming A (E) = Ay = 
= constant. Important results for this case have been obtained by WALLER [13], 
Wick [15] and others. The problem was also treated by the present author 
in a previous paper’, of which this paper is a sequel. The papers quoted give 
a fairly complete picture of the slow neutron density in the constant / case, 
including analytic expressions for this function at all distances r. 

When the variation of A4(#) must be taken into account, the difficulty of 
calculating the distribution function increases considerably in the general case. 
A lot of work has therefore been done instead on calculations of spatial moments 
of the distribution function [8], which is easier than the computation of the 
function itself. Although this gives some idea of the space distribution, an 


1 See Water [12], [13] or Marswax [8]. 


2 See Marsuax [8], for instance. 
° In the following this paper will be referred to as I, and formulas in it for example as I (1). 


15 209 


G. HOLTE, On the space energy distribution of slowed-down neutrons. II 


explicit expression of course would give much better information. But there 
exist only a few results beyond the age-theory formula for the distribution 
function. VERDE and Wick [14] have made such a calculation, but only for 
a very particular form of the function A(#). For the case when A (Z) decreases 
with decreasing H, Wick [15] has obtained some important results and derived 
the asymptotic behaviour of the neutron density for 7 > co. 

It is the purpose of this paper to extend the knowledge of the slow neutron 
density in the case of a variable mean free path. As far as possible we intend 
to calculate analytic formulas for the distribution function, thus completing 
the results obtained for constant A. The assumptions made are the same as 
in I, and it seems unnecessary to repeat them here. The starting-point as 
in I will be the Fourier-Laplace transform of the transport-equation, and 
developments and analytic continuation in the variables of this transform will 
also be used. Previous results for the case of variable A mentioned above, in 
particular those by Wick [15] will also be helpful. 

There is an essential additional difficulty connected with the use of the 
Fourier-Laplace transformation in the case of non-constant A. In fact, this 
does not seem to be useful at all unless special assumptions are made on the 
function A(H#). On the other hand, we have no close theoretical formulas for 
A(E£) (except for hydrogen) so that in any case this function must be taken 
from the neutron cross section measurements [1], [4]. By taking for A(E) a 
polynomial of sufficiently high degree, for instance, any accuracy desired may 
be obtained. It is still better, implying no restriction in practice to make the 
assumption first used by Marsuak [8] 


d, 
(1) 1(B)= 3, (7) Hone 


Here A, and d, are arbitrary constants, but d,=0. (Hy) = energy of the neutron 
source.) Further we shall always take dy) =0. The number of terms in (1) 
thus depends on the accuracy wanted. The case of constant A is obtained in 
(1) by taking only the first term, 2 =A ) and is a good approximation for 
sufficiently low values of Hy)(H < Ey). In many cases of practical interest a A 
variation with only two terms, those in Ag and A,, is a good approximation and 
means a considerable improvement over the assumption 2 = 4). However, no 
special assumption is made now on the number of terms in (1), although the 
actual numerical applications to carbon and water were performed with this 
simplified expression for / (1). 

If (1) is used for A(), the Fourier-Laplace transformation is applicable. The 
system of equations which results from the transport equation after the trans- 
formation and an expansion in spherical harmonics was first given by MarsuaK [8] 
and has also been used earlier [5] by the present author: 


(2) | (2 + 1) hn (7) fn (@, 4) + - dM (1) tats (Os. 4G) ae 


+ nfn—1(@, 7 + dy)] = ¥(@, 4) Ono; (n= 0.1; 2, ae 


The form chosen here for these equations assumes for the mean free path of 
the slow neutrons Ay that Ay #0, which is always fulfilled in practice. hy (y) 
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are the functions given by I (4), (5) and (6). With the special assumptions 
made about the neutron source q(w, 7) is given by I (3). From (2) the func- 
tions fy (@, 7), f:(@,7),... have to be calculated. Of these functions only fg 
is needed if we are merely interested in the total number of neutrons per unit 
volume, unit energy range No(r,#). This is related to the collision density 
Fy (r, #) which is ordinarily used, by an equation No (r, HE) = A(E) v1 Fy (r, E) 
(absorption is neglected, v = neutron velocity). The other functions /,, fo,... 
are only needed if we also want to study the distribution of the neutrons on 
the different directions of movement, but this will not be done here. 

Having calculated /)(@, 7), there remains for the computation of Fo (r, EL) 
the evaluation of the integral I (1), or 


1 mice +00 
(3) Fol, E)=—saag. | ane { dood? folo, 7). 
mt co 


Here x= log E,/E is the logarithmic energy variable and r is the distance 
from the source. For abbreviation we shall sometimes also use, as in I, the 
variable R= r/o, i.e. here the distance measured with the slow neutron mean 
free path as the unit of length. 

Before we start the calculations, we shall briefly sum up the main points 
in this paper. In Section 2 the physically important solution of (2) is determined. 
Here (4) is the essential equation. It implies that fp) can be written as a 
product of a function gy and the continued fraction I (2) from the constant 
4 case. This makes it possible to use the same method as in [ for one inte- 
gration where the essential contribution to the integral comes from the same 
main pole and then apply the saddle-point method for the other integration. 
Section 2 then contains a detailed calculation of the function gy by means of 
a power series in w?. The crucial equations used in deriving gp are (5) and (6); 
it is also shown how gp may be given as a closed expression, eq. (7) ... (16). 
The rest of Section 2 is devoted to a discussion of go. 

Section 3 contains a calculation of /y without special assumptions on / (EF) 
(except (1)). It uses essentially the same methods of computation as in I, 
and gives an expression for Fy where age-theory appears as the first approx- 
imation. Therefore it is useful only for sufficiently small distances and slowly 
varying A(£). 

The main part of the paper is Section 4, where Fy is calculated without 
the restrictions necessary for the expression given in Section 3. This is possible 
by assuming instead that A(#) is decreasing with decreasing H, which is the 
most important case in practice. Then the results by Wick [15] concerning 
the behaviour of fy can be used. They are first summarized, and then it is 
shown how they can be applied in combination with the general method of 
developments and analytic continuation in the variables, which is here, as in I, 
the main tool for computing the integrand of the saddle-point integral. Then 
F, is calculated, given by the main formula (37). This is completed for large 
distances by (48). Further we have for Fy formula (44), valid within an inter- 
mediate domain, but not applicable for a too rapid variation of A(#). Finally 
the complete asymptotic formula for r— co is obtained (49). A discussion of 
F, concludes Section 4. 
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Section 5 is devoted to applications. First some generalities about the nu- 
merical calculations, performed only in a special case of (1), are given. Then 
numerical results for carbon and water are presented, including also a com- 
parison with experiments for water. 


2. Determination of the function f, (w, 7) 
Solution of (2) 


(2) is an infinite system of difference equations in two variables 7 and n. 
For 4 = do, i.e. Ay =A, =::: = 0, 4 appears only as a parameter and we have 
only a system of difference equations in n. The solution of this system [12], [13] 
may be given as certain continued fractions as in I (2) for fy (, 7). In the fol- 
lowing this function will be called /% (w, 7) = ¢/No(@, 4) and in general Nn (a, 7) 
is the continued fraction hn(y) +-:- contained in I (2). Now it is easily seen 
that the function fp(w,7) has the same poles (among others) as the function 
{%) (w, y). If we put 


(4) fo (a, n) ie {i (o, n) Jo (a, n); 


this also appears from a method of calculating fp (@, 7) directly in the form (4), 
where gy(w, 7) 1s regular and #0 in particular when Re(y) = 1 (the real part 
of 7) and @ is near the main pole w (6) (9 =y—1) of f%(, 7). This 
possible way of calculation is a generalization of a method used by the author 
in an earlier paper [5]. We may always obtain approximate solutions of (2) 
by breaking off the equations at a certain index k, and putting /, = 0 for 
n> k, as the condition lim /f, (w, 7) = 0 must be fulfilled for solutions of physical 
n—> co 
significance. The last equation (n = k) then gives directly an approximate value 
of fz, which contains as a factor the first approximation of the continued frac- 
tion 1/N;. Putting this value into the equation (2) for n = k—1, we obtain 
a difference equation in 7 for the determination of f,—-1 in terms of f;-2. This 
equation is of a slightly more general form than the equation treated in the 
paper quoted ([5], eq. (15) and (16)) but may be solved using the same method. 
The result is the product of the second approximation of the continued frac- 
tion 1/N;z-1 and a certain development, essentially a power series in w?, broken 
off for consistency with the term in w?. This value of /;—1 is then put into 
the equation for n=k—2 in (2) and the same procedure is repeated until 
we finally get fp as the product of the (k + 1)th approximation of /%) and a 
development, determined up to terms in w?*. It is important to notice that 
this development contains the function hy only with the arguments 7 + dy, 
Hn +d, 4 +d, + dy, etc. This method of solution is always valid for suffi- 
ciently small values of |w| and in passing to the limit k— co we obtain 
the form (4) for /y, where go is obtained as a power series in w? which in 
particular does not contain hg (yn) (but ho(y + dy) etc.). This expression for fo 
is convenient for the evaluation of (3) because we have the same main pole 
@® = @)(9) as in I and may perform one of the integrations by the method of 
residues, obtaining an integral similar to I (17). This contains a rapidly varying 
exponential function, making the application of the saddle-point method possible. 
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The method above for the solution of (2) thus proves the validity” of the 
expression (4) for fy and some important properties of gj. Unfortunately, the 
method is not convenient for an actual calculation of fy except for the first 
few approximations. It does not seem possible to obtain the general term in 
the power series of gg in this way. But instead we may try first to derive 
the development in w? of all the function fy and then sum up the factor /%, 
thus obtaining gy. Such a development of fp has, in fact, been used earlier [8] 
for the calculation of spatial moments of the neutron density, but has not, so 
far, seemed to be useful for the computation of the neutron density itself. 
However, combined with the transformation (4) and using analytic continua- 
tion, this development is very powerful, as will appear in this paper. 

We now proceed to the solution of (2) with power series in w. We try 


n! 
(2 =a 1)!! k=n, n+2,... 


(5) fn (@, 9) = (— 1)" (1)* fnx (7). 


The factor in (5) is different from that used before [8], but with the choice 
made here the system of equations which follows from (2) and (5) has its 
simplest form. An easy computation shows that these equations are 


1 Ay hy 
(6) fala) = pan | ress. 2 fnss,e-a rth) +, frassnaly +d) + gro do 
hal) |b aco, v=0 Ao 
with 
n2 
ee | (ka, ner 2.0 4,8, 2 = 0, 1, 2, 8): 


This is a recurrence relation for determination of fnx(7), starting with the 
value foo (7) = @/ho(y). In particular, we may calculate as many coefficients 
in the development of f9(@,%), fo2(7), foa(7),... a8 we want, though with 
rapidly increasing difficulty. 

Of course, it would be desirable if we could replace successive calculation 
by a closed expression for the solution of (6). By introduction of certain 
operators and convenient abbreviations for occurring sums, this is, in fact, also 
possible, and it seems worth while reproducing this solution here. For that 
purpose we first introduce an operator P,, essentially performing a summation, 
by defining 


(7) Be Dale ay ee 


D (yn + d,), 


where ® is an arbitrary function of 7. This abbreviation is immediately 
suggested by the system (6) which now takes the simple form 


(8) ink (7) = Tn+1 Pr fa+1,k—-1 (7) a Py fn—1,4-1 (7) F q Ono Oxo; 
(k=n, n+ 2,...,%=0,1,...) 


The solution of (8) will contain multiple summations, obtained by a repeated 
use of operators P, on functions 1/hm(7). Therefore it is convenient to in- 


213 


G. HOLTE, On the space energy distribution of slowed-down neutrons. IT 


troduce the following functions of 7, depending further on an arbitrary number 
of indices 2, 7, k,... 


, Lily 
9 Dijr... sad con 
0 


We notice the simple relation, which follows at once from (7) and (9) 
(10) Pr Diju...(Q) = Dnise... (7). 


Now fnz(y) obviously is built up by sums of terms (9) multiplied with 
factors 7», with different combinations of values of the indices 71,7,k... The 
problem is to find the law according to which these combinations are made and 
then obtain an expression for the sum in a closed form. This can be accomplished 
by introducing operators Q, (v= 1, 2,...), operating on the row of indices 
aj9k... in the following way. For every operator Q, we insert in the row of 
indices a sign, say ;, which divides the row into two groups. All indices to 
the left of the sign ; are left unchanged by Q,, which only operates on the 
indices to the right of the sign. These signs always occur in the row of in- 
dices in falling order v from the left to the right, so that the sign j is furthest 
to the right. Now we define the effect of Q,: it performs a subtraction of 2 
from the y first indices to the right of the sign ;; further it transfers all 
signs ;,».1,...i7 one step to the right, for example 


QO; Disjuice = Dijzovet.. 200 Qe Digit jim = Dipak 21 ae ete 


We may then define a product operator Q,...Q.Q>, which of course means 
a successive application of the operators Qo, Q.,...@,, and in particular a 
power Q (r = 0) is defined in this way (Q; = 1, i.e. it leaves the row of indices 
unchanged). Sums and differences of such products, also containing ordinary 
numbers as factors, are then immediately defined in the usual way. 

Observing further that the operators P, and Q, commute, and putting 
k=n+2y, it is easily shown that the following expression gives the solu- 
tion of (8) 


n+1 My +1 ms+1 Mo+1 
(11) fn, n+2» (7) aa i| > mn, | SS Won Do | - lms | SY Tiny Qpeia—m| : 
- t 1 


m,=1 WO ail Mz=1 m=1 
v y—1 2 


. Om3+1—m: m,+l—m, 4 n+1—m, 
we | = Q"4 Q, D (q)a;n4+1, 21 ...mte—lintentend elon 


W v—1 v 
(ayy =" 0,1 2 one 


In evaluating this expression we have to take the parenthesis » with the 
summation in m, first, obtaining then a sum of n+ 1 terms of the same 
general type as (11), though according to the effect of Qn tlm with different 
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indices in D; only the term with m,=n+ 1 has the same indices as in (11). 
Then in each of these terms the summation in m,—; has to be performed. The 
procedure is repeated until the last operator Q, has been used and the summa- 
tion in m, has been made. Now the term in (11) with m,=n+ 1 is in fact 
the result of applymg the operator rn41Pn on fn41,n+1+2(—1)(y) if (11) is used 
for this term, and the other terms with m,=1...n are easily seen to be 
the result of Pn operating on fn—1,n-142»(y), when (11) is applied to this func- 
tion. This proves that (11) is the solution of (8), since in addition the initial 
‘aps for n = vy = 0 are fulfilled by (11) which then gives fo (yn) = ¢Do(n) = 
= g/ho (). 

Here we are merely interested in (11) in the special case n = 0, and it was 
only a step on our way to obtain the form (4) for fo (@,7) and the general 
expression for gg(@, 7). This may now be written as in (5) 


Ms 
(sn 


(12) 9o (@, 0) = )?” Joay (7). 


The ea ents of the development of /)(w, 7), obtained from (11) by putting 
A, = 42 = --- =0, are called /%)(n). The coefficients go2,(7) then have to be 


calculated successively from a relation 


1 y 
(13) Jo2v (7) = q | (q ) foe (1 as ho (7) Pi to) 1) 902(»—i) | z 


(yaa l Rare rs) (Goo (Nt) oe Lhe 


At first sight this does not look encouraging, but it is simpler than one might 
think and an actual computation using (13) is not necessary. In fact, fos» (7) 
and f%) (y) contain according to (11) and (9) the factor 1/ho(y). But even in 
ho (7) foo» (nN) (vy >0) there are terms with a factor 1/ho (4), and in hg (ny) f% (n) 
all terms (for » > 0) contain this factor, which is directly seen from I (2). On 
the other hand, it was proved in connexion with the discussion of (4) that 
Jo (w,) does not contain ho(7y) at all. Therefore in (13) all functions of the 
right member except the first one have the factor 1/ho (y ). As the remainder 
after the subtraction of these functions from hg (7) fog» (7) 18 a function without 


this factor, g2»(y) is obtained from the function oy foz» (y) by cancelling 


precisely all terms with the factor 1/hy (7). Therefore we only have to make 
some slight changes in (11) to have the general formula for gog, (7). All func- 
tions Dojx...(y) are multiplied by ho(n)/qg and in the summation in (9) all 
terms containing a factor 1/h (yn) are excluded. More completely, we introduce 


An; Ang 
a ve hee 
4 Bjx... (4) = pees’ 
qd ) Whore (7) 24, hj (q zi An,;) hy ( ++ dn, siz dn) eee 
ny,=0 
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where >” means that terms with 1/ho(y) have to be excluded, or, observing 
the first number in the row of indices which is zero, jk...rsOt..., the — 
terms corresponding to nj = m =--- = Mp = Ns = My = 0. Then we have from (11) | 


2 m,—i+1 ma+1 | 
(15) 9p, (M) = 14 | > fox | ae EY Oe | Sei gum a | 


M,—1=1 M,9=1 mM, = 
a Tt Bed i ee 1 1 


v 


. ents gene) Laie Neary eyetee (v= 1, 2, 593 


v—2 y—1 


We have, for example 


(15’) 9oo(7) = 1, Go2(") = 71 E10 (7), 


Goa (n) = 11 [72 i210 (N) + 71 ror (7), ete. (7, = 1/3, rg = 4/15, ...). 


In general, the main part of the expression for a function for a certain » 
can be calculated from the expression for »—1 according to the following 
formula, which is useful in actual calculations. For that purpose we first de- 
fine as a generalization of g.,(7) a function gy, ;, ,,(7) which is obtained 
formally by putting the indices 7j7...*78 before the indices of all the terms 
Eiu...(n) in the expression calculated for g)., (7). It is then easily seen that 


(16) Doo, (N) = 11 Jox(v—1),10 (2) + 72 9o2(>—1),12 (n) + 
é ade m, _9+1—m 4—m™, 
+ 112173 | > Tm, 3 | > Tn, 42-2 OP | Orta 
m,_3=1 m,—4=1 
v—3 y—4 y—4 y—3 


E (y)128, 5...9-Vi9~..2105 (v = 273 am 


This simplifies the actual calculation of g)., (y) considerably, though even here 
the operators QY occur. But in (16) the number of successive operations is two 
less than in (15). For » = 2 the last term vanishes, and we get the expression 
in (15’); for higher values we have, for example 


(16’) 9o6 (7) = 71 904,10 (7) + 72 Yoa,12 (4) + 117273 Ey23010 (4), - - - 


Some analytic properties of fp (@, 2) © 


Nothing has been said, so far, of the range of validity of the development 
(12). Taking Re{m} larger than the largest of the numbers 1 — d, (vy = 1,2...), 
it is easily proved that (12) is convergent for sufficiently small values of ||. 
By means of analytic continuation go(w,7) then may be calculated in larger 
domains of the y-plane and the w-plane as long as no singularities are met 
with. 

In the constant A case we have for /%)(w, 7), except for the poles, only the 
fixed singularities at w = +7 in the w-plane. Taking first the case Re {7 }—>0, 
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remembering that we then have hn(y)—>1, it is seen directly from (9) that 
the variable A case is obtained from the constant A case by substituting 


A(Eo) = 4m = YA, for do, giving then go(w, 7) singularities at w = + i , as 
appears immediately from (5), (9) and (11). But go(@,7) is singular at these 


points and regular in the interval between them {if “* <1) even for finite 
m 


y-values. This may be shown from the original equations (2) as has been done 
by Wick [15] by moving step by step to the left in the 7-plane, expressing 
/n(@,) by means of the values of these functions at the points 7 + dz, 
OS: Si PE Se CN CD SR OS Rare | 

For all values of 6 = 7—1 in the right halfplane of 0, avoiding in the o- 
plane only the points @ = + w, (6), + w, (9 + dk), +o, (0 + dy + di), ... (k,l... = 
=1,2...), (and mw = +7), where +.@,(0) denote the different poles of 
1%) (w, 8), go (@, 0) is seen to be regular in this way, except for the fixed 


singularities at @ = +7 ce The points avoided are the poles of fo (@, 0), of which 


a 

all except w= = @ (6) are poles of go(w, 0). For this function w = + wp (9), 
in particular, is a regular point as well as the point 6 = 0 in the 6-plane. 
Therefore, for the 6§-values quoted, the expression (15) for gg(w, 0) may be 
analytically continued in all directions of the w-plane as long as these poles and 
ho 
dm 

The position of the singularity nearest w = 0 determines the radius of con- 
vergence of (12). Now the values of importance here are real values of 0, 
making all poles arising from + wp, lie between =—z1 and =+ 2, and purely 
imaginary values of w in the interval between the fixed singularities. These 
are therefore of importance only for Am =>», which also is obvious for phys- 
ical reasons. Now let us consider the frequently occurring special case when 
all 2, >0. Then for 6 real =O and purely imaginary w, the power series of 
9o(@, 6) contains only positive terms (9).,(0) > 0). The singularities of the 
series therefore are situated at the points of intersection of the circle of con- 
vergence with the imaginary axis, and must be either a pole singularity be- 
longing to + wo (9 + dn) (dm = the smallest of the numbers d,, dz, ...) or the fixed 
fly 
Jeg 
and the A variation are very small, taking a fixed real value for 6=0, the 
poles will lie outside the circle of radius dp/~Am, and then the radius of con- 


the points +72, +7 are avoided. 


singularities -+ 7 depending on which is nearest w = 0. Except when dn 


vergence of (12) is exactly ‘a This is always valid, if @ is large enough. 
my 


For the purpose of integration the primary development (12) is sufficient in 
this case, as it is useful for all values of the variables needed. 

In the general case with arbitrary values of A, it is difficult to predict any- 
thing about the circle of convergence of (12) for real 0 = 0, as the coefficients 
may have arbitrary signs, and it probably has to be estimated numerically in 
each special case. Possibly we then have to perform an analytic continuation 
to the point iAg/Am (if Ag/Am <1). 
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3. Calculation of the neutron density at small and intermediate distances 
when 2(E) varies slowly with E 


When 4(£) is a slowly varying function the evaluation of (3) may be done 
in quite the same way as in I, as long as r is not too large. The result is 
a function Fy (r, #) containing corrections to the age-approximation, which, so 
far, is the only general expression known for the neutron density for small r. 

As in I, we may perform one of the integrations by the method of residues. 
If we integrate w in this way, we have to consider all the poles mentioned 
before and to sum the corresponding residues. But in restricting ourselves to 
the slow neutron density, it is sufficient here, as in I, to take into account 
only the main pole w =a )(6)=7B(6). It is true that we now have many 
new poles, contained in gg (w, 9), but these also are all easily seen to be “‘sec- 
ondary”. In fact, even the contribution from the pole which is the most crit- 
ical one w= @ )(0+dm) is seen from (2) to be only of the magnitude 
exp (—dma) as compared to the main pole term. If we suppose x2 10, for 
example, this is sufficient for neglecting it, as dm in general is of the magni- 
tude = 1. 

Calculating now the residue of fy) (@, 6) at w = @ (8), and writing it, as in 
I, for small |6| as a power series in 0, this series is obtained in the follow- 
ing way. It contains first the factor resulting from / determined by I (18), 
(14) and (16). In this paper we shall use for I (14) the notation A)(6), and 
we write oo for the coefficients of its power series in 6 which are calculated 
by I (16). Then we have the factor resulting from gj), which is obtained by 
development of g.,(8) around 6 = 0, and inserting in (12) w =72B(6), which 
gives the function Gy (6) 


(17) Go (8) = go [@o (8), 0] = 1+ 9,0 + 9207 +--- 


Here 9;, 92... are numerical constants, but contain all A, as parameters. The 
main numerical difficulties in solving the problem with 4 variable lie in this 
computation of g, and they increase rapidly with increasing values of », and 
of course also with increasing number of parameters 4, which are used in the 
formula (1) for A(#). Concerning the function Gy (@) we know that it must 
have a singularity for the value of 0 giving B(0@) = Ao/Am. The root of this 
equation (for Ag/Am<1) is called 0m. On the other hand, this must be the 
first singularity on the real positive 0-axis, as the substitution w = w,» (0) in 
go (@, 8) makes all the singularities deriving from the poles of go (@, 0) move 
farther to the right. Therefore G (9?) may be calculated in the interval 
0=0<O6m on the real axis, either directly by the series (17) or by making 
an analytic continuation of the series in addition. In the special case A, > 0, 
if the A variation is not too small, we have the radius of convergence 0, for 
(17) directly. 
We now define 


(18) A (9) = Ao(4) Go (9), 
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given for small |@| by a development 
(19) A(0)=1+ 040+ %0?+---, 


just as in I (14), and for other values of 0 as the analytic continuation of 
this series. The formula which results for Fo(r,#) is then I (17) or I (36) 
(after division with A>) and the remaining integration in 6 (or ~) may be per- 
formed by the saddle-point method as in I. As soon as dAg/Am< 1, however, 
there is an important restriction in the applicability of the formulas of I, 


dependent on the new singularities at = +i * of fy (w, 6) in this case (or, 


Am 
physically, dependent on the importance of the primary neutrons when their 
mean free path is larger than that of the slow neutrons). In fact, I (17) and 
I (36) are now valid only if the contribution to the integration in w from the 


integral around the point wm = +7—* along a cut in the imaginary axis (sim- 
mm 


ilar to the cut from wm =+ 7 in I) can be neglected. 

But this is possible only when the real values of 6 of importance in the 
evaluation of I (17), ie. the saddle-point 0;, is considerably smaller than 0m. 
This also appears directly from the fact that A(6;) in this case is singular for 


6; = 0m. When 6; is near 6m the singularity at w = + i gives the essential 
m 


contribution to the integral and the evaluation of (3) is then possible only with 
a close knowledge of the nature of this singularity. Now if r and Am/Ao are 
not large 6; may be taken as in I as the root of the equation « — R B’ (6) = 0, 
and therefore R/x = r/Aju <1/B’(Om) is the condition for I (17), I (36) 
and the formulas derived from them to be useful. The extension of this do- 
main depends essentially on the value of Am/Ay. If this quantity is not too 
small (say 22), the condition is roughly r<22/¢/O? Am, obtained by tak- 
ing the first term of I (11). In many cases, slowing down in water, for ex- 
ample, this is a very small distance, and other methods of calculation must 
be used. 

When ,r is sufficiently below the limit above, we may take over directly 
the main formula I (40), divided by 4, for Fo(7,#), and also I (41), (42), 
(43) and (44). Now for very small 7 we must of course expect these formulas 
somehow to contain the age-approximation. They also give the expected Gaus- 
sian behaviour, but in the same form as in I, exp (— r?0?/4/¢ x), whereas the 
age-approximation for the variable 2 case contains instead of A) a mean value 
A, I (49), or in this case with the A-variation (1) and for sufficiently large x 
given by 
3 A 


9 tices sens aia oe | 
(20) iz = ia = (ae + 2) = il. Pan FN | 


(k,l =0,1,2,...,4=1=0 excluded), 


formally corresponding to an increase of x to vw =a+ a. 
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In order to solve this apparent inconsistency, we must first examine the 
approximations made in age-theory. This theory is better the larger M is, and 
in fact may be considered as an asymptotic theory for M— oo. In this limit 
the age-theory formula must follow from the integral (3). Following now the 
evaluation of (3) to the next step, I (17), and making M-— oo, this solves 
the problem. For in this limit Ap (0) = 1, corresponding to the constant A case, 
but the new factor G)(6) of A(@) in the variable A case, has non-vanishing 
coefficients in the limit M-—- co. This depends on the limiting form of hy (6), 
ho (0) >20/M for Moco (but hn(6)>1 for n=1), making the term 
r Eyoio10...(8), which is the largest one in gp, (4), behave like M’, whereas 
wo (9) +—60/M, so that g.,(0)A0. Indeed, Go(6) has the limiting form 
exp (a0) for M->oo, which follows because the age-approximation then is 
exact, remembering (20) and I (17). Of course this can also be proved by a 
direct calculation of G (0), taking into account the facts mentioned above, 
but the proof is not very simple and the reader will be spared it. 

Now we can make use of this property of G)(#) even for finite M. It is 
natural to separate the factor exp (a0) from G) (0), writing 


(21) Gp (0) = 62° G0 (8), A(0) = e7° A’ (8). 
The new functions may be determined for small |6| by the series 
(22) Go(0)=1+ 910+ 9207 +---, A’(O)=14+ 0404+ 2 6?+---, 


and these coefficients are obtained from g, and «,, computed before (17), (19), 
in developing exp (a6). Then we can still use I (17), only substituting x’ for 
x and A’(0@) for A (6), and similarly in all the expressions evaluated for Fy (7, £). 
Now, this is advantageous in every respect. Firstly the age-theory formula 
now appears as the first approximation of our expression for Fy (7, Z), and 
the next approximations contain in a natural way the corrections to age-theory. 
On the whole, we obtain an extension of the domain of applicability of our 
formulas in this way. For the coefficients «, are smaller than a, at least for 
y small, which is helpful in developing Fp (7, EZ). The increase of x to w’ =a +a 
(for a>0, which is always valid in practice) is also advisable in this re- 
spect, as 1/a’ is essentially the parameter of these developments. This increase 
is also important because it decreases the saddle-point 0;, which extends the 
formula to larger values of r. For 6s is now determined by the equation 
av’ — RB’ (0s) = 0 (for r and &Am/Ag not too large), similar to I (33), and the 
maximum value of + by the condition R/a’ = r/Aga’ <1/B' (Om). 

We may now take over all the formulas I (17) and I (33), (34)... (44), 
changing x to xv and A(@) to A’(0). It seems advisable to reproduce here: the 
simplified form of I (40) for small R/z’, as in this case, when « in general 
is not small or decreasing with increasing » (in contrast to the constant 4 case), 
the developments for /(r,#) as given in the earlier papers mentioned [5], 
[6] are best written in a somewhat different way. We obtain, either by a 
partial development in I (40), or directly from I (36) 
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(23) Fo (r, #) = 


Q- 3h, usieeel I On OF 3.08 301Y 
deh jr (a tyh” AV REE hoy it tagegeet 3 | a 


alee 4 


U! () [2-24 BAe .| aoa i ey ...|— 


(a’ t)'2 De rele. 24 


Su (8 eere || Ce 


Here 6, and gs = + V6, are the saddlepoints in the 6- and g-planes, respec- 
tively, given together with t by I (37) (1#>2’,R=r/d). The exponent in 
(23) then has the form 


R? C7 (1 + 24) 
3a’ (1+ 7)? 


Faay) us. 2’ 0, — RB (6) = —| + Rd(p), 


d(y) =—O(B.y + Bay’ + Bag? +>). 


Further we have in (23) 


(25) Uy =~A (P)=9+ uP + ag? +---, 
and 
(n) 
(26) Sn fle Mel in Soe weet oe Ge 


(a'r)? on} (a’ rt)" 

In (23) the term —R6(qgs) in the exponent, as well as the terms dn and y, 
are only small correction terms, negligible for sufficiently small R/a’. Then 
we also have t ~ 1, and (24) transforms into the Gaussian exponent — R? C?/4 a’. 
In this small domain around the source (roughly for R/«’ <1/CV38,), (23) 
may be simplified (see earlier paper [5], formula (31), taking «—> a’ and 
% — a). Now in (23) we must always have 6;<Om, implying an upper 
limit for R/z’, but also, for the validity of I (37), roughly R/a’ <2V2/0V3 py. 
Of these two conditions the first one is the more restrictive when Am/Ag is not 
near unity. Then (23) is the only simplified form of I (40) which is applic- 
able to the variable A case. For a very small A variation, (23) is completed 
at larger distances, compatible with the condition 0; < 60m, by the approx- 
imate I (44), with 6; taken from I (38) (making the usual changes). The 
condition for all these formulas to be useful is, as usual, that the first terms 
of the developments decrease sufficiently rapidly. 

The corrections to age-theory appear in (23) partly by the deviation of (24) 
from the value — R?C0?/4.2’, and partly by the deviation of the development 
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| | | Fie) 9 DE aaa 

in (23) with the main term na U (qs) ei Gs 32 (2' t) from the 
value RO/22’. These deviations both increase rapidly with increasing R/z’. 
The first one is of the same kind as in the constant A case and is also the 
more important for a slow variation of A. The second one also appears in 
the constant A case though it is then of little importance. In this case, how- 
ever, it is of dominating significance as soon as the A variation is not very 
small. For a large variation of 4, in particular for small M, the transforma- 
tion (21) is not sufficiently effective. Then the development of (23) or I (40) 
may be useless even for 7 = 0, which indicates that age-theory is a very poor 
approximation in this case. 

These formulas for the neutron density, giving the corrections to age-theory, 
have been derived for the case of a slow variation of A. They are applicable 
only up to a certain distance 7, determined by the condition 6; < Om, de- 
pending physically on the increasing importance at larger distances of the 
source neutrons with the mean free path Am >Apg.. If Am is not much larger 
than A’), implying a sufficiently low value of Hy, and when /A) and M are not 
small, these formulas are sufficient for practical purpose. Otherwise, it may 
be necessary to complete them. This will be done next, though we are re- 
stricted then to the practically most important case of a mean free path that | 
decreases with decreasing energy of the neutrons. 


4, Calculation of the neutron density when 7 (E) decreases with decreasing E 


New information about fp (@, 0) 


The extension of the formulas for the neutron density to larger distances, 
or to the case of a rapidly varying mean free path, requires knowledge of the 
kind of singularity possessed by fp (@, 0) at @ = + %Ao/Am (for Ap < Am), (or of 
Go(9) at 0 = 6m). This seems to be difficult to obtain in the general case. 
However, in the case where A(E#) decreases with decreasing #, Wick in an 
important paper [15], has succeeded in calculating the behaviour of fy (@, 6) 
at this point, and has thereby derived the asymptotic dependence of F% (r, E) 
on r and w for rocco. Using here this information about fp (@, 6), and com- 
bining it with our previous knowledge of this function and also making use 
of analytic continuation, we will now calculate the neutron density at distances 
outside those already attained. In the special case roo, the formula then 
obtained automatically shows the asymptotic dependence predicted by Wrox. 

Now, in this case we get from the integration in in (3) an essential con- 
tribution from the cut along the imaginary axis from the point w = + idg/An- 
Therefore it is convenient to take instead the opposite order of integration, 
calculating first the integral in 6 by the method of residues. The poles of 
{%(w, 0) then are the inverse functions of w= «,(0), say 0=0,(m), i. 

= 0,[w,(9)], in particular for the main pole 0 = 0)(w). Further we make a 
shght change of variables in (3), putting 


(27) @ = tk, 
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and denote the main pole henceforth as 4)(k). This is obtained by taking the 
inverse function of k = B(@) or k® = [B(6)]?, and can be calculated for small 
|&| as a development in k?, from I (11) 


(28) Oo (k) = A/C? +b, h (C4 +---, 
and for larger & from I (27) 


EN 1 k—hy ha (k—ho)? ; | 
(29) Oo (k) [5 + m i a 


When the values of the coefficients in I (11) and I (27) are known numer- 
ically, the coefficients of (28) and (29) may also be rapidly computed numer- 
ically. These developments are also of use in the previous calculations for the 
value of Om (k = dg/Am). We also point out here the simple approximate for- 
mula for 6 )(k) contained in I (32), valid for 6 real, positive and not too 
small. 

The result of the integration in @ is then easily obtained, neglecting as be- 
fore the contribution from all “secondary” poles. Denoting simply by A (4) 
the function resulting from (18), (19), after substitution of (28) for 6, where 
A (k) is given for small |/| as a development 


(30) A(k)=1lt+aP+a,k*+---, 
we have 


8x7 Er 2 


(31) Fy(r, Z) = fa e~kR+260(%) 4 (k) 09 (k). 
Ly, 


This formula is similar to I (17) and I (36), and may also be evaluated by 
the saddle-point method. The path of integration Lz, which originally is the 
imaginary axis, may be transformed to pass through a saddle-point k; along 
a path of steepest descent. When A(k) can be considered slowly variable, 
ks may be taken as the root of the equation «0o(k) — Rk = 0, corresponding 
to «—RB’'(@)=0 before, having the value k; = B(6;). With the simple ap- 
proximate values B,(6) of B(6), I (31) and I (382) ks is equal to B, (6s), 
giving 

2RO*{2 47) ba 
(32) ks 321 + oP and k,=a lee 


The first equation is obtained from I (31), (37), and is valid for small R/a; 
the second is obtained from I (32), (38) and holds for larger R/z. 

But these values of ks are useful only for sufficiently small r and for A/dm 
near unity. Then (31) may be evaluated in quite the same way as I (17) or 
I (36) and gives equivalent developments for Fy (r, #). The same transformation 
as in (21) may naturally also be used to extend the range of validity to larger 
distances. However, we are going to use (31) only when the previous formulas 
for Fy(r, Z) are insufficient, and the transformation (21) must be replaced by 
a more powerful one. 
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Now Wick has shown, that for 0 — co to the right in the 0-plane, fy behaves 
like (1 — kAm/Ap)?-? for k > Ag/Am, where @ is a constant <2. Here @ may be 
determined as the lowest eigenvalue of a certain Schrédinger equation (WICK 
[15], eq. (105), cf. also Appendix. in the present paper), of which the higher 
eigenvalues give additional terms to the term quoted of the same type as this 
one, although these are of vanishing importance at k = Ay/Am. Then we can 
move step by step to the left in the 6-plane, solving (2) successively for de- 
creasing values of 6, as was described before. The result is that as long as 
the poles are avoided, /) behaves in the same way for kA /Am even for 
finite 6-values. 

This singularity of course belongs to go (ik, 6), which shows this behaviour 
in particular for real 620. Then, however, we may have a pole in the k- 
plane between the origin and the fixed singularity if 6 is small. But after 
inserting 0 = 09(k) in go(ik, 6), all these poles move beyond k = Ag/Am, and 
the resulting function go (ik, 9(k)) is regular in the interval between k = 0 
and this point. With the variation of A(#) now assumed, in any case if all 
4, > 0, it was proved before that (12) has the radius of convergence Ag/Am 
when 6 is real and chosen sufficiently to the right in the 6-plane. This is 
valid also after the substitution 6 = 0)(k), if the variation of A is not too 
small. But this function g(7k, 09(k)), given by a development in k?, multiplied 
with a similar development of Ao (9% (k)), resulting from the residue of the 
factor {%, is precisely A(k) (30). Thus A(k) is regular in particular for 
0<k<AglAn. Further the development (30) of A(k) in general converges 
within the whole of this interval, if Ag/Am is not too small. For &— do/Am, 
A (k) behaves like (1 —kAm/A9)~*, taking « = 2—o(a>0). We therefore put 


L (k) P (k) 
(1— kdm[Ay)*  [L— Ue Am Ag)? 1"” 


(33) 2 A(k) = 


where the second expression is taken because A (k) is an even function. 

L(k) and P(k) are obtained as power series in k. We know that for 
0 <k<Ag/An in particular, L(k) (and P(k)) is regular, and that it is continu- 
ous also at k=Ag/Am. We cannot assume regularity at this point too, but 
only that L(k) consists of a regular part + terms singular like (1 — k Am/A )2i~2, 
containing the other eigenvalues 0;(i = 1,2...). However, these terms all 
vanish at k= A/dm. Then to the first approximation, the behaviour of (33) 
around this point is determimed by knowledge of the value L(Ag/Am) = Lm. 
This value may be obtained from the power series in k for L(k) = A (k) 
(1—kAm|Ao)* (or for P(k)) and then either by putting k = Ap/Am directly in 
this series or by first using analytic continuation. For this purpose various 
methods may be devised, and an account of this will be given later in con- 
nexion with the reproduction of the numerical applications. The essential fact 
here is that such a calculation always is allowed because L (4) is continuous at 
Ag/Am and is possible in practice by computation of a number of coefficients 
@,42,... of A(k). Now, it seems in general that for the next lowest eigen- 


value 0; we have ej3—@>1. Then the behaviour of A (hk) at the singularity 


may be obtained to the next approximation, too, calculating graphically, for 


* Similar transformations (with a= 1) can also be used for the study of 6° (k) (28) at the 
singularity k= 1 (cf. I). 
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example, the value Lm = L'(do/Am) and developing the regular part of L (kh). 
The graphical study of L(k) also gives a rough estimation of the coefficient 
Im of (1 —kAm/Ap)-2 in the complete development of L(k) around k = Ag/Am. 


Calculation of Fo (1, E) 
Making the substitution (33) in (31), and putting 
(34) K(k) = L(k) 60(), 1 (k) = P (k) 60 (hi), 


we now proceed with the evaluation of (31). When 7 is not too large, we 
take the singularity at A%»/Am into account best by writing it as a logarithm 
in the exponential (cf. BAzRwatp [2]). Thus we have for the integrand in (31) 


(35) T(k)e%™, with W(k) =—kR + 205 (k) — a log [1 — (kAm/Ao)?I.- 


This transformation is similar to (21) but more effective because it sums up 
the main singular part of A(k). We have used here the second transforma- 
tion (33), because P(k) is an even function (7'(k) odd) and easy to handle 
numerically. Its power series converges satisfactorily at & = Aj/Am and may be 
used directly without any preliminary treatment. 

We may now use a saddle-point integration in (31). But instead of (32) 
which is restricted to small r or Am/Ag, we take ks as the root of the equa- 
tion (at least approximately) 
iy , 2 oldie Ae )e 
(36) YW" (k) R+ x 0 (k) + 1— (bdml Ay)? 


0. 


Ly, is chosen as a straight line through k,; parallel to the imaginary k-axis. The 
resulting expression for F(r, #) is obtained and may be written essentially in 
the same way as the general formula I (40) 


Q 3h, 
(2V2)? Br E 


(37) Fy (r, BE) = he eM [irae] eee 


Here the development is exactly the same as in I (40), and is therefore not 
written out. The quantities d, «, mand n; in this development now are 
abbreviations of 


(38) d =[P" (k,)/2]-"2, & = dW" (hy), 
_ gi (he) 


PO (hs) 
Mi. = = 
a! 


CEO per Ge Bae 


a! 


M; is obtained by I (41). If the saddle-point ks is determined exactly, which 
then in general has to be done numerically, we have ¢=0 and M; = m. 
Otherwise we may take 


bo 
bo 
OU 
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RO 
ie 2 x" 


(39) ks > a =at+a(Cdm/d)?, 


as long as this makes e small (i.e. for sufficiently small 7). This gives for 
Fy(r, Z) a Gaussian main term for small 7, but with a coefficient of —r? in 
the exponent different from that of the age-approximation (x > 2’). 

(37) is our main formula for Fo(7,#). It is applicable for much larger 
values of 7 and Am/Ay than (23). As in I (40) e ought to be small, and m, 
M; and » are regarded as small quantities of the ith order of magnitude 
(~ 1/(Vz)'). The development in (37) as in I (40) is essentially an asymptotic 
series in 1/Vz, containing only even powers, which are calculated up to the 
fourth order. The formula thus contains three successive approximations of 
Fy(r, Z). The zeroth approximation is obtained by taking only the term Mo 
in the development, the second approximation by adding 


— My (15 n?/16 + 372./4) —3M,n,/4—M,/2 to this term (cf. I (40)). 


The complete development calculated in (37) and I (40) gives the fourth ap- 
proximation. Then we can get an account of the accuracy of (37) (as of I 
(40) before) by studying the values of the successive approximations. This has 
to be done numerically in each special case. Then it appears that the simple 
zeroth approximation already gives comparatively good values. The second 
approximation in general changes the value by about 10—50 %, whereas the 
fourth approximation in general implies a change of a few % only. Although 
the fourth approximation may be needed for great accuracy, the second ap- 
proximation therefore ordinarily is sufficient in practice. Due to a consider- 
able cancellation of terms, this also holds when the order of magnitude rela- 
tions for each n; and m; are no longer strictly valid, which is the case for 
increasing r and Am/Ay. In fact, a check of (37) with numerical integration 
of (31) along Ly shows an excellent agreement even in such cases. At the 
same time, however, (39) no longer makes ¢ small and is therefore useless, and 
we have no simple general relation ks; = ks(R, x, Am/Ao), although it may be pos- 
sible to calculate for given values of x and Am/Ay. Thus in a small inter- 
mediate domain we cannot obtain the general analytic dependence of F on R, x 
and AmldAg, though (37) of course still is useful by giving first certain values of 
ks within this domain and determining then the corresponding values of 7 
from (36). For still larger r, when ks > Am/Ap we could again derive a simple 
approximate analytic expression for ks by developing Y (k) around k = Am/Ag. 
However, although (37) even for 7» co shows the correct dependence on 7 
and « (cf. formula (49)), the development then seems to become uncertain 
and it is convenient to modify the method of integration as will be described 
below. We may also expect (37) to begin to fail in accuracy for increasing 
values of Am/”~o, roughly when (CAm/do)? is of the magnitude =x. Possibly 
numerical integration of (31) along Ly; then will be needed for more accurate 
values of Fo. 

For large values of r we are faced with the problem of a saddle-point inte- 
gration when the saddle-point approaches a singularity of the integrand. Then 
it 1s convenient to develop around the point of singularity instead of around 
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the saddle-pomt. In fact, problems of a similar kind have been treated before, 
by Barrwatp [2] for example, though in the case of a single pole. Here the 
singularity is a branch-point, too, which necessitates some changes of the 
method. We take in (31) 


(40) k = dglAm + 9, 
and develop 4 (k) 
(41) Oy (k) = Om + yp Om + y? Oml2+ ++, 69 = 6 (Ay/Am). 


According to (36), the condition for ks to be near dg/Am then can be approx- 
imately written R—2z Om > «Am/dg. Then we use the first transformation (33), 
multiply according to (34) with the derivative of (41), taking the following 
approximate expressions of L(k) and K (k) 


(42) L (k) = Lm + Lin y + ln (—p Al Ag)®-°, 
K (k) = Ko + Kip + Koy? + ky (— ye. 


Here Lm... are the constants mentioned before, the new constants Ky... are 
calculated from them by (34) and (41). Kg, is not completely determined in 
this way, but still we can take it as a very rough estimation. Now (31) takes 
the form 


” 


x 6, 
Q-3hy (40\" 6 a —(R=2 9) 0+ pees 
(43) Fo(r, E) 872 2 Er ae ev"m inp : dpe - 


Ly 


“(—y)"* (Ko + Kip + Kay® + ky (— p22). 


I, is the path of integration in the y-plane, cut along the real axis from 0 
to co. L, may be chosen for convenience so as to make the integrand in 
(43) decrease as fast as possible from w= 0, only subject to the restriction 
that it always lie to the left of the cut. This choice will depend essentially 
on which of the two main terms of the exponent, —(R—2 On) yp or px O0m/2 
is the dominating one. The first one decreases more rapidly for p increasing 
along the real axis (since R—xOm>0), the second one for y leaving the 
origin along the imaginary axis, as we always have 0;,>0. Therefore two 
cases can be distinguished 


1) wAm/Ay < R—2 Om 5 V0 Om/2. 
This requires a sufficiently large value of 0m, say 0m210. Therefore this 


case applies only to moderate values of Am/Ag. We take L, as a straight line 
to the left of, but infinitesimally close to the imaginary y-axis. The terms in 


y, p>, p',... of the exponent are then expanded, and the result is a develop- 
ment containing integrals of the type 1/7 | dy &” (— p)-* (B = « Om/2). These 
i 
v 
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integrals are connected with J-functions in a very simple way. This is seen 
in detail by the transformation y =— 2’, determining the radical by a cut 
through the negative real z-axis. Then the integral above transforms into a 
well known contour-integral for the J™function, and has the value 


We obtain the following expression for Fy 


(44) Fo (r, £) = 
Q- 3h, e%m—*4m (* a 1 [*, LB lel: =| 
ib SS a = se 
8A he Er ed Hy Bert V2 67/21. 39m Ko 
a i be Qo+3 Fe nes Ky By, Qo+3 ss K, al ae 
EE 1 fe test gf 
“ Q1— Q 0 
(a Om/2) 
Here we have 
Oy pe), pte) ey ayy ae 
(45) Q_, (243) u (r=) Py (é) 1+ oh ! 4! 92 t > 
. —— 
2 2 
5 Bei 
Pr eel Mec ne SAC Re e = (R—29m)|Vx 6/2, (R=1]Ap). 


see noe oe 


A condition for case 1) is e <1, but, in fact, (44) is applicable also for slightly 
larger e, if only the terms of the development in 1/Vx6%,/2 in (44) decrease 
sufficiently rapidly. This main condition for the formula to be useful is ful- 
filled for «>1 and 0, > 1, the latter of which restricts the applicability of 
(44) to moderate values of Am/Ag. 


> a Anlas 
2) (R—2£60m) 4, > Vx 0%,/2. 
> Ky/Ko 


In this case the path of steepest descent from y= 0 lies near the real po- 
sitive p-axis for small |y|. For increasing values of Re{ yw}, however, this 
path leaves the real y-axis and passes to infinity roughly parallel to the im- 
agimary axis. Therefore we can still take for one part Ls of Ly two straight 
lines parallel to the imaginary axis, starting from points infinitesimally close 
to a point k= ky in the interval Ag/Am< ks <1, one on each side of the cut. 
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A convenient choice for ks is the value given by (32), for |then the exponential 
part of the integrand decreases as fast as possible from k = ky along L,. The 
other part of L,, L,, joins these points with each other. L, is a contour along 
the cut, infinitesimally close to it, which encloses the origin in the p-plane 
(k = Ag/Am). Thus we have in (43) is = (i So if . In this sum the last term, 
Bee lace he 
the ““branch-point” contribution, always dominates the first term, because it 
contains the main term e~"/4m, whereas { behaves like e~® = e~" for 7 co, 
Ls 
For the calculation of f the terms in y*, y®,... are expanded, and then we 
ie 
: ; 1 , 
obtain a sum of integrals of the type ; fay e-?¥ (— y)-* (vy = R—a# Om). Now 
Le 
such an integral, where the path of integration is extended to y = + co in 
both directions, is precisely the contour integral for the J-function encountered 
above. From this value we must here subtract the contributions from the 


paths between the point k=; and the real, positive infinity. Then the inte- 
gral quoted has the value 22 (y)*-!v(s, u), where 


(46) (0) = Fe|! reset 


Here Q(1—s,u) is an incomplete J™-function 
(47) Q(1—s,u) =f dtett*, with w= (ki—Ag/Am) (R—2 Om). 
Using these notations, (43) takes the form 


5 Om — 1/4, ’ 
Q 3h e (R209)? Ky | 2 — 0, ¥) — 


Reteeis BA da oes. Br 


fat oh i] 2 a 
Ky (R— x Gin) @ (R— x On)? 
For (2 a, &) | pth 
+ K,(R—2 6,2 + + Fos (r, E). 


Here Fos(r, #) is the contribution to the integral from Ls, and may be cal- 
culated essentially as a saddle-point integral (if ks is not too close to Ag/Am). 
It contains the same main term as I (40), exp (—Rky + @0(ks)) and has 
the character of a correction term in (48). It will be negligible when w is 
sufficiently large. At the same time the incomplete J/™-function (47), of magni- 
tude exp (—vu), can be dropped in (46). Actually it seems possible in most 
cases to use (37) or (44) and (45) up to a distance where w attains a value 
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say ~ 5. If we use (48) for larger distances only, we can substitute 1/I’(s) 
for v(s,«) in (48) and drop the term Fos(r, £), which simplifies the formula 
considerably. Possibly there will be a small intermediate domain where neither 
(37) and (44) nor (48) will give accurate values of Fo(r, #). Then we can 
interpolate the neutron density at these distances with sufficient accuracy for 
practical purpose. 


The asymptotic formula. Discussion of Fo (1, E) 


We notice that for very large distances (48) is simplified, giving the follow- 
Ing expression 
Q Om — tm 3h, hn Ko ; 
bad, EG | eet eo) 


(49) Fo (r, #) 


(49) is the complete formula for the asymptotic slow neutron density. It shows 
the dependence of Fy(r, #) on x and r predicted by Wick [15]. Further it con- 
tains the constant yw, depending on M and on the rapidity of the A-variation, 
and which can be determined numerically as has already been described. But 
from (48) we also obtain direct information on how large r-values we have to 
take for the asymptotic formula to be valid. The condition is r/Ap > x Om, 
completed for large Am/Ay with R > K,/Ky, which essentially implies r > Am 
(a similar, but weaker condition must be fulfilled also for (48) to be valid). 
For r/A) < % 0m we found before that the roughly Gaussian behaviour of (37) 
was dominating. For increasing 7, when r/Aj9 — & Om > G@Am/Ao, it appears from 
(44) and (48) that we have already the same main term exp (% 0m — 1/Am) in the 
expression for the neutron density as we have in the asymptotic formula, 
though multiplied with a different factor. But this shows anyhow that even 
at these distances the rapid source neutrons of mean free path Am are of domi- 
nating importance for the slow neutron density, even if the complete asymp- 
totic formula is not valid before the limit Ay x Om (and Am) has been consider- 
ably surpassed. 

In a discussion of the neutron density the value of Ag xv Om is a character- 
istic quantity. In general, 6 (as well as 0m) decreases with increasing values 
of Am/’o (M = constant), but increases with increasing values of M (Am/dAy = 
= constant). This gives the general result, to be expected of course for phys- 
ical reasons too, that large values of wx, 49, M, and small values of Am/A 
(near to 1) all are factors which extend the validity of the roughly Gaussian 
behaviour and make the domain of the asymptotic formula and the related 
formulas (48) and (44) apply to larger values of 7. If Am/Ag is not too small 
(roughly = 2), this appears directly from the approximate value of Ag x 0m ~ 
= 2225/AmC*, for large M further simplified to Max/3/3 Am (I (45)). As was 
mentioned before, this is quite a small distance for water, for example, which 
means that (44), (48) and even (49) have great practical importance in this 
case. But for carbon, for example, this quantity is much larger, which then 
makes formula (37), (23) or I (44) most useful, though even here (44) and 
(48) are applicable at distances of practical interest. These facts will appear 
in detail later in connexion with the reproduction of the numerical results. 
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A discussion of the neutron density as complete as in I is more difficult in 
this case, but the most important facts have been given above. For small 
r-values and a slow variation of A, the discussion in I on the whole is ap- 
plicable (e > 2’). But for larger values of r and Am/Aj, there are also striking 
differences between the constant A case and the case treated here. If we take 
a constant value of x and study thé space distribution, this may be described 
as in I as a quasi-exponential decay, though now with the decay-length de- 
creasing towards 2». But this behaviour is now in general more marked and 
begins earlier than in the constant A case. As it has been pointed out before 
[15], the energy distribution of neutrons at a certain distance 7 has quite an- 
other form here than in I for large r-values, as it tends to a limit independent 
of r for r—>co. This also implies that the energy distribution function is not 
so rapidly decreasing with increasing values of HK as in I. Its dependence on 
M and on An/Ap is determined essentially by 0 and the variation of Om was 
quoted above. Thus for large Am/’yp we have roughly 0m ~ (Ao/AmC)?. 

Further, we have in (44), (48) and (49) the constants 9 and mw (Ko). For a 
roughly twofold variation of A and a moderate M-value o@ is near unity. @ in- 
creases for increasing Am/Ag and varies only very slowly with M. Finally 
/ Increases with increasing M, roughly as M for large M, and with increasing 
AmlAg; it tends to a limit independent of Am/A, for large values of this quan- 
tity. These facts appear in detail from the numerical results. 


3. Applications 


Methods of numerical calculations 


A numerical calculation of the first corrections to the age-approximation, by 
means of a computation of g, and g,, for example, can always be made with- 
out undue labour, even if several terms are used to describe the function A (EF) 
(1). On the other hand, it is natural to restrict the number of terms in (1) 
so as not to be unduly large, giving an unnecessarily refined formula for A (#). 
We have to take into consideration the uncertainties of the experimental meas- 
urements of A(#), for instance. Further, as we are now concerned with larger 
values of Hy than in the constant A case, other general assumptions made, like 
considering only S-scattering between the neutrons and the nuclei, are not ful- 
filled with the same accuracy. Anyhow, we must expect our formulas for the 
neutron density to have less precision now than for those small values of Ep 
for which I was applicable, and therefore it is not worth while taking too 
many terms in (1). In this paper we have used for the numerical calculations 
merely a two-term variation of A(#), capable of describing a linear, and also 
a more general kind of monotonic 4 variation 


(1’) A(E) = Aq + Ay (;,)- 


In many cases of practical interest this is sufficiently accurate. 
The numerical calculations at first aim at a computation of the coefficients 
91,92,--- Of Gy(@) in (17) or ay, a,,... in (19). Besides, the knowledge of 
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B(0) and 69 (k) is essential; but this has been obtained before in I, as well as 
the factor A»(8) in (18). Go(@) is the result of the substitution w = 71B(0) im 
the series (12) for go(w, 0) (9 =7—1) and a development of the coefficients 
Joov (9) around @=0 (y= 1). The main numerical difficulties lhe in this cal- 
culation of gog,(6). Now, in the case treated here, the summation in (14) 


contains only terms with nj, nz,-:- =O°or 1. Any function Ejx...(0) with » 
indices jk... occurring in go(w, 6) can therefore be written as a polynomial 
of degree v in A,/Ao 

(50) | Bye,..(6) = > (ial2a) HD. ._(0) 


Here the “components” EH (6) of Hyx...(6) (« = 1, 2...) depend only on 6. 
In the future the value of this component, given as a development in 6 around 
the origin will be denoted by Ei °) . What we finally want is thus the devel- 
opments EY >) of the pomponeate of all the functions F;;..,(6) occurrimg in 


the expressions evaluated for go2,(9), (15). These coefficients may also be 
written as in (50) 
2» 


(51) GJo2r (8 (Ay [Ao)* y of; (9). 


2 
When the component es (0) («= 1,2...2¥) is developed around the origin 
the resulting power series in @ is ae ge). 

For the basic functions h; (0) in all these calculations the numerical method 
of obtaining developments around 6=0 was given in I, Appendix B (for 
M>1). As these could be calculated rather rapidly up to a high degree in 9 
they can be used even to get the developments around the points 6 = ud of 
h; (8), though of course only up to a lower power in 9: This has actually been 
done in the case of carbon. For M = 1 the original formula for h;(@), which 
is simple in this case, had to be used, I (4) and (5). Finally we invert these 
series, calculating the developments of 1/h; (0) around 6 = wd (u = 0, 1, 2, ...), 
and denote this by 1/h” (i and mw not both = 0). 

If more than the first or the two first terms in (17) are wanted, it seems 
best to use the following recurrence relation for the calculaticns, obtained 
directly from (14). 


(52) LG ae ss ra Eye u + en Bie Hebe mle, w= 0, 1, 2, . +) 
with the initial conditions BO” = 1/hM, HL = 1/h“+, defining = = 0. 


Here a generalization EN “ of the function Ey °) | defined before, hag been 
introduced. This is Wi Colin of (Ay/Ao)* in ‘the expression (14), when 
this coefficient is developed in a power series in @ around the point 0 = ud 
(7=1+ 4d), (u=0,1,2,...), with the summation in ,, mg... extended 
(for ~>0) to all combinations of the indices of summation. Thus for the 
function Be i given as a power series in 0, the restriction concerning the 


summation in (14) has to be used only for ~=0. But in computing all the 
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functions H\%,) which enter into the functions gog, (0), successively by means 
of (52), eis restriction is taken into account automatically by the formal 
requirement 1/h = 0. 

Now, when the basic developments in 6 1/h™ (2, 4=0,1,...) have been 
calculated, we also know the initial quantities H%” (x = 0, 1, =), Vo. 
for the application of (52) (because the last index in all functions nae (0) i mn 
9o2»(9) always is zero). We compute first H%” and then we have to calculate 


E&Y, Be!) (x =0...3,4=0,1,...). We could go on in this way until all 
the sae pepater a Bee wanted have been obtained, and then insert them into 


the expressions evaluated from (15). But instead we define developments 
g3/°, obtained by substituting Ee me f6r 1 ae everywhere in the expression 
for af ©; and similarly we define functions Gig ten ee Then the calculations 
may be done more rapidly by making use of (16). In fact, for the formation 
of these functions or, on the whole, any function which is a linear expression 
in the quantities H%”) a formula, completely analogous to (52), holds. g%”, 
for example, (15’), is obtained by applying this formula to the development of 
11 (72 EY) + 7, Eo), with i= 1. For the computation of g%”), (16’), we then 


use this formula to get r, 9%) and reg and apply it once again to the sum 


of these quantities and calculate 7, Gee t) + r299¢4- Then we only have to add 


11273 E4,, which must be computed separately, using (52). For higher func- 
tions gi" (v= 4,5...) (16) is used and the main part of such a function is 


rather rapidly computed in the same way as above from the values of the 
immediately preceding term. The rest of formula (16), the operator expression, 
which remains to be evaluated, causes a little more trouble, but even here 
several terms can be put together to form expressions g/%". ;.., derived from 
lower values of ». 

The numerical difficulties are considerably greater here than in the case of 
constant A, and it does not seem possible to compute as many coefficients 
a, (g) aS in I. But fortunately, this is not necessary, even if we want the 
same accuracy as in I, except when the A variation is so small that we are 
near the limiting case 2 =) = constant. For when A, >0 and A,/dy is not 
too small, the development (30), which has to be analytically continued to 
ke = AglAm = Ag/(Ag + 41), already has this point lying on the circle of con- 
vergence, which is advantageous compared with the cases of analytic continua- 
tion met with in I. But the number of coefficients may be reduced also’ 
because the same accuracy as in I, for example, is now useless for reasons 
mentioned above. In fact, the actual calculations, including the cases of carbon 
and water, were performed only up to g;, but this seemed to be sufficient for 
practical purposes. Then we had to start with the comiputation of the values of the 
functions h;(9) and their first p derivatives at the points 0 = ud for all values of 
1, p and «(u~=O0) compatible with the equations w+7+2p<10 andi+ p<5 
(in general: <2n, and <n, respectively, if n coefficients g, are wanted). In this 
way the basic quantities 1/h\ were obtained as power series in 6 up to 6”. 


EeY (4=0, 1,2) is then calculated for w=0, 1... 8, BGs and £uY 
(x =0,...3) for w=0,1...7 etc. Hach step in these computations consists 
of the formation of the sum of two products, each containing two factors in 
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the form of power series (52). The coefficients of the resulting power series 
are then sums of a number of products and are rapidly obtained numerically 
by the use of an electric computing machine.t This computation is quite 
analogous to the numerical evaluations of equations I (21), ... (25) and is per- 
formed in practice by the same methods. The numerical calculations which 
remain for obtaining g,...g5 and a... a5, when all g% («=1,2... 2», 
y=1,2...5) finally are known, are also of the same type. 

In a same way as in (51) a, can be written as a polynomial of degree 2y 
in Ay/Ao 


(58) OY oral yf dif 


The values of the constants a,, are the final results of the first and the largest 
part of the. numerical computations. Of course explicit formulas of the first 
few of these constants, containing the values of h;(6) and its derivatives at 
6=0,d,2d,..., can be derived. This was done in an earlier paper [5] for a, 
and ad . It is also possible to calculate for a, developments in 1/M, valid for 
large M, similar to I (45), (46). Thus we have, for example 


2 5d if 10d 
2 ahh OU eee a mie) sete carl, oa fe 


A IbdiMe ia ade nail ade 15 M Pea en 
a uae ee one le 
re Mts 15M 
But these developments are not at all as useful as those for C, 61, 62, ... and 
019, Gg9,--- given in I, because the convergence is much poorer and is still 


worse for increasing values of » and x. (54) shows the expected asymptotic 


2 
behaviour for Moo, remembering that we now have a = 2A, ae 1 (2) : 
diy 2d\Ay 

(20), (21). 


Proceeding now with the calculations, we obtain from B(6) the coefficients 
of (28) for 6)(k%) and insert them in (19). Instead of a, we then have coeffi- 
cients a, (30), also containing A,/4) as parameter, but which could also be 
written as in (53). Now, if A,/A4) is not small, both a, and a, will have large 
values, increasing rapidly with », and therefore it is convenient for the nu- 
merical computations to take a new variable ¢ in (30), writing O (¢) for A (k). 


2v 
(55) k= tdg/Am = tag/(Ap + Ax), C(t) =1Ltatt+ettt-:-:, oy = ay (3 
m 


* Facit, mod. ESA, Atvidaberg, Sweden. 
* These formulas may easily be derived also for the general case (1). 
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Then the series for C(t) has the radius of convergence = 1, at least when A, /Ao 
is not small. In the actual calculations this was confirmed by the facts that the c, 
were close to unity, and the power series for C(t) is of the same type as the 
simple series 1 + (2+ ¢4+- 

For the calculation of K,, K, ..., We now first compute the eigenvalue 
o(=2—a) using the variation method with a linear variation function (Ap- 
pendix). We could then, according to (33) consider the function DL (tAg/Am) = 
= C(t) (1—t)*, which is continuous also at t=1, and try to calculate this 
function and its derivative at this point. But the power series in ¢ for this 
function has alternating signs with the coefficients near unity, of the type 
1—t+f—.---. It is true that such a series can be summed for ¢t = 1 if it 
is first conveniently transformed, using “Euler’s transformation”? for example.+ 
This was also actually performed in some cases. But there are simpler methods 
available for this computation. The reason why the development of C (¢) (1 — t)* 
has the inexpedient form quoted, is that this function is still singular and in- 
fiite at ¢ = —1, as C(t) is an even function of ¢. Therefore it is much better 
instead to study the development in ¢ of the inverted function (1 — t)~7/C (d), 
which is finite and continuous for —1<t<1. This series ought to be of the 
same type as 1/(l—t+2#...)=1+4#, and it was actually found that the 
coefficient of ¢ was near + 1, whereas the other coefficients were smaller and 
decreasing. The convergence of this series seemed to be satisfactory for ¢ = + 1, 
so that it could be used for determination of the functional values at these 
points. Then we also have a certain control of the convergence of the series 
from the sum at ¢ = —1, which must be zero. 

This way of calculating LZ, and Ky seemed to be the best one, but a third, 
and still simpler method was also used, and the values resulting from the three 
methods were compared. This makes use of the development in ¢ for the 
function P(tdg/Am) = C(t)(1—#)*, which is contmuous for —1St< +1. 
This is rapidly computed, as only even powers of ¢ occur, and the series showed 
satisfactory convergence at t= +1 (it ought to be of the type (1+ @+ 
+ ¢#*+ ---)(1—#) =1, which also was confirmed). The sum at t = 1 divided 
by 2¢ gives Lm (and Ky). The graphical study of L(t4p/Am) in the neighbour- 
hood of ¢=1 then gives a rough value of L;, (and K,), and finally in this 
way an estimation can be made of Im (and ko). 


Application to carbon 


This case was treated in I for A constant, and for the values of B(0) and 
Ay (6) we refer to Tables 3 and 4 in this paper. According to the measure- 
ments of the neutron cross section for carbon [1], [4], the treatment in I covers 
all values of Hy roughly up to 0.1 Mev. For higher Hy, roughly up to 2.5 Mev, 
the following expression gives a very good approximation of A (£) 


(56) A(E)en = 2.65 + 2.35 E (Mev), thus d= 1, dy = 2.65, A, = 2.35 Ey (Mev). 


1 See Harpy, “Divergent series’, Oxford 1949, (1, 3, 3), and (1, 3, 4), page 7, observing 
that in the latter formula 6, = 4" a), if a, =(—1)"a,. (A a = Ar-1aqi — An-1q)) 
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Then we have taken the carbon density as 1.60 g/cm*; for any other value, 
(56) of course has to be multiplied with a constant. First we give the expres- 
sions (28) and (29) of the inverse function 6) (k) of k = B(6), obtained from I, 


Oo (k) = 2.2214 k? + 1.4116 k* + 1.1556 k® + 1.0372 £8 + 0.9741 KO + os, 
(57) 4 0) (k) = (0.3 — 1.8500 u + 1.897 u? — 7.68 u? — 5.0 u4 — 52 u® — ---)—1, with 
u = k— 0.8326, 


the first of which is used for insertion in (19) and besides for determination 
of Om roughly for 2p/Am < 0.6, and the second used for 0m for Ag/Am = 0.6. Then, 
however, I (32) and I, Table 4 give a simple approximate formula which also 
may be used. For the derivatives of 6)(k) more accurate values can be ob- 
tained by differentiating the inverse function k = B(6), using all the coefficients 
of I, Table 3, when 0, has been determined before. 

The result of the main part of the numerical calculations, the constants 
Qx, 18 given in Table 1 (for ao see also I, Table 3). 


Table 1 
Oy x 

bk 

aN? 1 2 3 4 5 
0 — 0.30435 0.05258 | —0.00520 0.00024 | —0.0000 
1 2.2722 —0.0532 0.0269 —0.0023 —0.000 
2 0.6413 3.2809 0.8236 0.1476 0.012 
3 1.9873 4.4416 2.098 | 0.576 
4 0.3148 3.8588 6.221 3.801 
5 1.2570 | 6.379 8.930 
6 0.1441 | 3.057 | 9.946 
7 0.7087 | 6.037 
8 0.0646 | 2.069 
9 | | 0.881 

10 | 0.029 


Then, for every value of Ey within the interval 0.1—2.5 Mev (roughly) we can 
obtain the basic constants a,... a; by a very rapid computation according to 
(53) and (56). For low values of Ey and small 7, when (23) is applicable we 
can then first compute the constants a, by (21) and (22), with a = 2A,/A» + 
+ 3(A,/A)®. For the high values of Hy and for larger distances we transform 
(19) into (30), using (57), and after the calculation of @ (a) we apply (33). Then 
we can use (37), (44) or (48), and some constants occurring in these formulas 
are given in Table 2 for three different Hy. Rough values of these constants 
for other Hy may be obtained by interpolation, at least for e, Ky, Ky, ..., which 
are not easily computed directly. 

In Table 2 the last, or possibly the two last figures of the values for 0 have to 
be regarded as uncertain, as well as the last figure of Ky, w and K,. For Ke 
and ko only one figure is significant. For the two highest values of Ej we have, 
very roughly, o,—@ ~ 2. Thus the term in hp in (44) and (49) is approximately 
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Table 2 
Mera pn | ep fand md we 
| 
2 0.3155 1.916 7.39 0.7646 | 1.68 15 40 100 14.5 
1 0.7691 3.628 14.1 0.4615 | 1.54 10 30 50 9.3 
0.5 1.621 7.49 40.0 |—0.1519 | 0.44 Liz 


of the same importance as the term in K,. For the lowest value of Ey quoted 
an estimation shows that K, and Ky, are of the same magnitude as Ky. 
However, these values are omitted in the table because they are very rough; 
besides, they are of no practical importance in (44) or (49). As a demonstra- 
tion of the calculations of Ky,... we take the value Hj) = 2 Mev. Then C(t) 
(55) has the following development 


C(t) = 1 + 1.6583 @? + 2.178 ¢* + 2.57028 + 2.87048 + 3.106 2° + ---. 
Further we have 


1/L (tA/Am) = (1 — t)-2/C (t) = 1 + 1.2354 t — 0.2776 #2 — 0.560 # — 0.141 ¢ — 
— 0.112 4 — 0.016 #8 + 0.003 27 + 0.003 2 + 0.006 2 + 0.001 2° + 0.000 #1 + ---, 


which converges well for |¢| = 1, though a more accurate value of o (a) probably 
would make the convergence still more rapid. It seems to give a reliable value 
of Im (and Ky), and by a graphical study of the sum near ¢= + 1, even of 
Im (and K, and K,). The corresponding developments for the lower values of 
Ey) quoted show a somewhat slower convergence, but still seem useful. But as 
was pointed out above, for decreasing Hy accuracy of the values Ky, Ky,... 
becomes less important. 

Finally we study the validity of the formulas for Fo(r, #) in this case. For 
Ey =2 (23) is not even applicable for r = 0, but instead (37) is useful up to 
fairly large distances. (48) is valid only when 7 has surpassed the value 5.082 cm 
sufficiently, but such a distance may still be of practical interest. The asymp- 
totic formula (49) already begins to be approximately valid for r > 30a cm, 
roughly. For the lower values of Hy (23) may be used for small 7, though (37) 
has a still larger range of applicability. The characteristic limit Ay Oa of r is 
larger now, 9.612 and 19.82 cm for Hj) = 1 and 0.5 Mev, respectively, which 
makes (44), (48) and (49) still less practically important in these cases. The 
asymptotic formula (49), for example, requires distances roughly > 60a and 
120 z cm, respectively, to be approximately valid. 


Application to water 


The slowing-down of neutrons in water has been studied theoretically by 
several authors. Fiiiccr [3] has given a simplified treatment of this problem, 
which, however, is based on the diffusion equation and therefore must be con- 
sidered as insufficient. VERDE and Wick [14] have developed a numerical 
method of calculation which is apt only for this case because it is based on 
the fact that A) is small for water, which is taken into account by putting 
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Ag = 0 in (1’) (€= 0.5). But this approximation is very rough when the slow 
neutron density is wanted, even if it is possible to apply corrections afterwards. 
The application of the general treatment given here, with both terms in (1’) 
taken into account from the start, seems more reliable. 

However, we have so far studied the slowing-down only in moderators con- 
taining one single element. Even if the oxygen atoms do not contribute very 
much to the scattering and slowing down in water, they ought to be taken 
into account, at least approximately. The simplest way of doing this is by 
assuming that the scattering cross sections of hydrogen and oxygen vary with 
E in the same manner, so that their ratio is constant and independent on #. 
Denoting as before the total mean free path in water as 1(#), and the mean 
free path in hydrogen, i.e. in water, when only the scattering of the hydrogen 
atoms is taken into account, as A4g(E), we have A(E)/Aq(#) = x = constant, 
and we only have to take everywhere in our formulas (cf. [13]) 


Cn (9) = %Cnx (yn) + (1 — x) eno (m), 


if cna (yn) and cno(y) are the functions cp» (7), I (5), (6) for hydrogen (M = 1) 
and oxygen (M = 16), respectively. In this sum the main contribution always 
comes from hydrogen, as x is near unity, and in particular for n > 0, because 
then ¢cno(y) becomes very small compared to cnx (7). 

Now, from the measurements of the cross sections [1], [4] we can take the 
following expression for the total mean free path (units for E and Ey Mev) 


(58)  A(E)om = 0.65 + 2.09 (E)", thus A) = 0.65, d = 3/4, a, = 2.09 (Ey). 


Except when # is near a resonance value for the cross section of oxygen, this 
formula is quite a good approximation in the interval 0—2.5 Mev, and can be 
used even for higher energy values, though with less accuracy. Making a com- 
promise between the values for small and large #, we can take, with good 
approximation in the interval quoted, x = 0.87. 

First we have to calculate the function B(@), using the methods developed 
in I. But since an approximation A= A) = constant is not useful even for as 
low a value of Hy as 0.01 Mev (giving roughly Am/A 9 = 1.1), this function is 
not needed for large values of 6, and the analytic continuation to @ > co with 
the determination of the constants of the asymptotic formula for constant A 
is of very little practical use. Therefore, only the constants 3h,, OC and Bn, 
I (11), An, I (27), and some values of a, b and c, I (32), were calculated. Thus 
we obtain 


(3hy = 1.244, 
(CO =—3 1.04950, py — 0.48408, Bo = = eOOie etc. Bs hems Bur 

(59)? 4 0.05819, 0.17818, — 0.4556, 0.7920, —1.218, 1.781, —2.55, 3.64, —5.2, 
respectively, and 


ho = 0.660, hy = —0.127, hy = 0.035, hg = —0.009, hy = 0.002, ... 


; 1 In these constants, as well as in other constants given, of course only two or three 
figures may be considered as physically correct. For mathematical reasons, however, more 
figures must be taken and treated as if they were correct. 
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The values h», refer to a development I (27) around the point y,; = 1.5, 6, = 2/3, 
derived from the first development by means of yp = 6, 0) = 1/6 as intermediate 
pomt. For the constants of I (32), used for the saddle-point 0, in I (also useful 
for Om), we can take 


(60) a=0.981, b= 0.360, c = 0.456, for 0.3565 2.0, or 
D502 h20.63 (U.6e Bias 17. bige). 
The inverse function 6) (k) of k = B(6) has the following developments 


*{ 9% (&) = 9.90789 k? + 0.7980 k* + 0.8270 k® + 0.8071 k® + 0.8261 k1° + ---, 


61 
eh 69 (k) = (1.5 — 7.89 u + 17.0u2? — 36 u? + 60ut + ---)-, wu = &—0.660, 
the first one used for the insertion in (19) and also for determination of Om, 
roughly for Ag/Am<0.55, the second one for Om for larger values of Ag/Am 
roughly up to 0.8, which seems to be sufficient in practice. The first one may 
also be used for the derivatives for small Ag/Am. These can also be calculated 
more accurately from B(6), differentiating implicitly, as more coefficients were 
computed in (59). 
Now we proceed to the constants a,,. They are given in Table 3. 


Table 3 
Oe 
nv 1 2 3 4 5 
SX 

0 | —0.88984 0.71730 | —0.4943 0.2442 0.016 
1 3.8888 —0.6937 0.879 —1.069 1.33 
2 1.1187 11.319 4.415 1.443 —1.84 
3 7.2830 33.157 25.91 7.27 
4 1.2175 31.455 101.76 100.6 
5 10.923 117.84 319.3 
6 | 1.3196 61.05 415.7 
7 15.00 278.2 
8 1.430 102.1 
9 19.62 

10 1.55 


For very small values of Hy the development in 6 with the constants obtained 
from this table is useful directly. But in most cases the transformation to the 
variable k by means of (61) is needed, and the formula (37) for Fo (7, £) is 
the most important one. Formula (48), however, is also of great practical 
interest in this case, much more than for carbon, for example. The same con- 
stants as for carbon above are given in Table 4 for four different values of Ko. 

For the accuracy of the values of the constants 9, Ky... we have the same 
limitations as for carbon. An estimation of the next lowest eigenvalue gives 
the very rough value 9, ~ 2.5 for all the values of Hy quoted. The develop- 
ments used for the determination of Ky, ... show good convergence, at least 
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Table 4 
Eo | OF, a On, : 0 Ko \ Ky Ko ko i 
2 0.0226 | 0.296 2.06 1.0001 | 0.207] 27 10 5 1.65 
1 0.0538 | 0.477 2.44 0.9212 | 0.306 | 2.8 8 4 1.67 
0.5 0.1197 | 0.781 3.38 0.7907 | 0.403 | 2.6 7 3 1.59 
sah) 0.309 1.60 as 0.4841 | 0.42 2 6 6 1.16 


for the three highest Hy). As must be expected in general, the lowest value 
of Ey gives the poorest convergence. Therefore for Ey = 0.2 the values of Ky, 
Ky, and ky are only very rough estimations. 

As was mentioned above, formula (37) for Fo(r, #) is the practically most 
important one in this case. In spite of the fact that water ought to be an un- 
favourable case for the application of our formulas on account of the particularly 
high values of Am/Ao, (37) seemed to be more useful than one might expect. This 
appeared from a numerical study of the successive approximations of (37) and a 
check with numerical integration of (31). Using (37) and (48) (and for very small 
Ey (44)), possibly combined with interpolation through a small intermediate domain, 
this seems sufficient for the calculation in practice of Fo(r, H) at least for the 
values of Ey) quoted. Except for small 2), the functions obtained show a roughly 
Gaussian. behaviour only within a very small domain around the source. The 
range of practical interest is characterized by a transition from the Gaussian 
to the exponential dependence on 7, a sort of quasi-exponential decay with a 
decreasing decay-length — A». The purely exponential decay predicted by the 
asymptotic formula (49) ought to begin to hold even at comparatively small 
distances. For Hy) = 2 Mev, for instance, a distance 7225 cm should be suffi- 
cient for this formula to be approximately valid for slow neutrons. Finally, a 
numerical comparison of our formulas with the age-approximation reveals that 
this approximation gives considerable errors even for 7 = 0 (values for Ey = 1 
of the magnitude ~ 20% too small, for example). But as this approximation 
predicts a slower decrease of F(r, #) for small r than our formulas, it will 
also give values of Fy) which are too large within a small domain near the 
source, before it finally begins to decrease much faster than the actual expo- 
nential decay. 

In conclusion, some comparisons with experiments will be given. In fact, 
several measurements of the density of slow neutrons in water have been made. 
Most of these concern thermal neutrons, but there are also many measurements. 
of resonance neutrons, which are the main interest here. Unfortunately, the 
most used source, an Ra—a-— Be source, has the energies of the primary neu- 
trons almost continuously distributed within an energy range 1-10 Mev (roughly), 
so that an application of our formulas requires a superposition of terms cor- 
responding to all the different values of Hy. This is rather laborious though 
possible in practice of course, if we know the primary spectrum of the Ra —a— Be 
source. But the knowledge of this seems to be incomplete. Recent measure- 
ments by Trucuer [11] indicate two maxima, at about 4 and 6 Mev, respec- 
tively, in this spectrum. Qualitatively we may say then, that on account of 
the values of Hy in question Fy (or rather 7 Fo) must be expected to show an 
exponential behaviour already at distances ~ 15-20 em from the source (condi- 


240 


ARKIV FOR FYSIK. Bd 3 nr 14 


0.3 


0.1 


ay 10 15 


Fig. 1. Hyp = 0.22 Mev. 
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tion for (48) to be valid). The decay-length ought to be roughly ~ the mean 
free path of the primary neutrons corresponding to the larger maximum Ky = 6 
Mev, or approximately ~ 9 cm, which also is indicated by experiments, by 
Rusu [10], for example. Further the decay-length ought to be slowly increasing 
with r (in contrast to the case of a mono-energetic source), on account of the 
increasing importance for large r of the higher energy components in the primary 
spectrum. 

Photo-neutron sources are in general better known and are often almost 
mono-energetic. For this case a comparison was made with measurements by 
O’Ngat [9]. In figures 1 and 2 the function r* Fy) is shown for D neutrons 
(= In resonance neutrons, # = 1.44 ev), calculated for Hy = 0.22 Mev and 
0.9 Mev respectively, using (37) and (48) (for Ey = 0.22 also (44)) and assuming 
the same value of Q in both cases. These are the primary energies assumed 
for photo-neutron sources of RdTh + D,O and RdTh + Be, respectively, and 
the measured neutron density from these sources, brought to the same maximum 
height as the theoretical curves, is also shown in the figures. As the source 
had an extension of roughly 4 cm, and r is the distance from its centre, a 
displacement of the measured curve from the theoretical curve to the right 
must be expected. This displacement ought to be near 2 cm in the neigh- 
bourhood of the source, and then decrease with increasing 7. Taking this into 
account, the agreement between calculated and measured curves may be con- 
sidered good. 


Appendix 


Calculation of the constant @ = 2—a. Wicx [15] has shown that @ can be 
calculated as the lowest eigenvalue of a certain differential equation ([15] eq. 
(105)), formally corresponding to the radial Schrodinger equation for a two 
dimensional particle. For the application to water this equation must first be 
generalized to a medium containing more than one slowing-down element, which 
is possible when the mean free path /;(#) in the component 7 varies with E 
in the same way as the total mean free path A(#). Then we take 4/A; = x = 
= constant (2 = 1, 2,..., > “i= 1). If we assume for simplicity only two 

v 


different nuclei (%; + %, = 1), with masses M, and My, neutron masses, respec- 
tively, the equation is 


(62) a : £ + iE Re 3) + Aye“ + Age “a — Fh oth T() =10, 
with 

(63) AS 2 — Bee Bey. ee 

Here we must always have y>0, given then by 

(64) y = Eh’ (Eo) /4 (Eo) = 2X4, Ay! Am 
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if (1) is used. For the A variation actually treated numerically, (1’) we obtain 


(64’) y = day lam: 


At first sight (62) looks suspect, but it is in fact symmetrical in M, and Mg, 
peda : ; M yeas 
which is easily seen by the transformation & = o|/ The generalisation 

1 


to more than two components is obvious. 

The lowest eigenvalue @ of (62) can be calculated by the variation method. 
Wick [15] has used the eigenfunction exp (—so?/2) with s as the variable 
parameter for a rough computation for carbon. Using the general equation (62), 
the minimum condition ‘is | 


47° st 8 8 


ay" te ae ene Crees eked 
(65) + Aya + we ae ee ine @ MJ, +s 


However, this method was found to give too bad values of 9. But (65) was 
still used for calculation of a rough value sg for s. In order to obtain a better 
value of o we then tried a linear variation function 


(66) U(c) = Smee, with z= oVs/2, 


n=0 


making use in this way of the approximate knowledge of the eigenfunction, 
which results from (65). In (66) sy is a constant, taken roughly as the root 
of (65), c, are independent variables, which have to be calculated from the 
usual system of linear, homogeneous equations. This also contains the eigen- 
value as a parameter, which is the root of a secular equation 


(67) Det. | Ann’ | = 0, (m, nm’ = 0, 1, are), 
where the elements of the determinant have the form 


(n +n’)! 


(68) Han’ ea “on+n' +2. EE gare Fh Ay pias +1 @+ 


Gener) 49" ; 
| : +n’ +1)I- 
So n+n' So Mi rena ) 


Here a, and a» have the form given in (65), inserting sg for s. @ is the para- 
meter to be determined by solution of (67), related to the lowest eigenvalue 
o of (62) by the equation 


M, 
PG ore : 
(69) o <1.5 2 Din 
Then ®,, is the largest root of (67). This must be determined numerically by 
inserting trial values of ® in (67). These values can be estimated with the 
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help of the roots of (67) for the determinants of the first and the second order, 
which can be computed directly. In practice it was only necessary to take 
two or three trial values of ®, and then calculate the root ®n by mterpola- 
tion. However, using a determinant of the sth order in (67), the value thus 


obtained is still only an approximation ®” of ®,,, always smaller than the 


. . . (2) 
correct value. Therefore, the successive approximations Dy’, On, ... were also 


calculated and the convergence studied. Then it seemed sufficient to use the 
sixth order determinant in (67), and the last figure in the value chosen for Om 
was obtained by estimation from the successive approximations. For the nu- 
merical evaluation of (67) for each trial value of ®, the method devised by 
JAMES and CooLIDGE [7] was used. In this way we also made an estimation 
of the second eigenvalue 0. 

For the application to carbon the formulas are simplified because we have 
%, = 1, x, =0. As an example, for Hy = 2 Mev, we took so = 0.7 and the suc- 
cessive approximations of @®, are: 0.251, 0.260, 0.308, 0.309, 0.3149, 0.3151, 
and the value 0.3160 was chosen for ®,. For water we took hydrogen as 
nucleus 1, thus M, = 1, x, = 0.87, and M, = 16, x, = 0.13. The terms corre- 
sponding to oxygen then only appear as small corrections in (62), (65) and 
(68). Taking then Hy) = 2 Mev, for example, we used sp = 2.6, obtaining for 
®» the following successive approximations: 2.47, 2.47, 2.52, 2.52, 2.525, 2.527, — 
and for ®m the value 2.531 was chosen. 


SUMMARY 


The stationary space energy distribution of slow neutrons in an infinite, 
homogeneous slowing-down medium, where the mean free path of the neutrons 
A(E£) depends on the energy H, has been studied. /A(#) is assumed to have the 
form 4(£) = >A, (E/E) (d, = 0, dy = 0), where Ep is the energy of the neutron 


source. This is supposed to be an isotropic point source. Firstly an analytic 
expression for the collision density Fo(r, #) has been calculated (r = distance 
from the source) which is valid for a slowly varying A(£) and at sufficiently 
small distances and contains the first corrections to age-theory. Secondly, the 
case of practical importance when A() decreases with decreasing # has been 
treated and general expressions for Fy (r, #), covering essentially all values of 
r, have been calculated. The term which dominates the behaviour of Fo(7, Z) 
undergoes a transition from a roughly Gaussian to an exponential dependence 
on 7. In particular, for roo the function calculated shows the dependence 
on r and £ predicted by Wick. The numerical constant in the asymptotic 
formula is also obtained, as well as detailed conditions for the validity of this 
formula. Numerical applications have been made to carbon and water. Com- 
parisons show good agreement with experiments for water. 


The problem of the slowing-down of neutrons was suggested to me by Pro- 
fessor Ivar Water, and I wish to express my gratitude to him for his en- 
couraging interest in my work and for valuable discussions and suggestions. 
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